home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dlasq3.f < prev    next >
Text File  |  1996-07-19  |  22KB  |  821 lines
  1.       SUBROUTINE DLASQ3( N, Q, E, QQ, EE, SUP, SIGMA, KEND, OFF, IPHASE,
  2.      $                   ICONV, EPS, TOL2, SMALL2 )
  3. *
  4. *  -- LAPACK routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       INTEGER            ICONV, IPHASE, KEND, N, OFF
  11.       DOUBLE PRECISION   EPS, SIGMA, SMALL2, SUP, TOL2
  12. *     ..
  13. *     .. Array Arguments ..
  14.       DOUBLE PRECISION   E( * ), EE( * ), Q( * ), QQ( * )
  15. *     ..
  16. *
  17. *     Purpose
  18. *     =======
  19. *
  20. *     DLASQ3 is the workhorse of the whole bidiagonal SVD algorithm.
  21. *     This can be described as the differential qd with shifts.
  22. *
  23. *     Arguments
  24. *     =========
  25. *
  26. *  N       (input/output) INTEGER
  27. *          On entry, N specifies the number of rows and columns
  28. *          in the matrix. N must be at least 3.
  29. *          On exit N is non-negative and less than the input value.
  30. *
  31. *  Q       (input/output) DOUBLE PRECISION array, dimension (N)
  32. *          Q array in ping (see IPHASE below)
  33. *
  34. *  E       (input/output) DOUBLE PRECISION array, dimension (N)
  35. *          E array in ping (see IPHASE below)
  36. *
  37. *  QQ      (input/output) DOUBLE PRECISION array, dimension (N)
  38. *          Q array in pong (see IPHASE below)
  39. *
  40. *  EE      (input/output) DOUBLE PRECISION array, dimension (N)
  41. *          E array in pong (see IPHASE below)
  42. *
  43. *  SUP     (input/output) DOUBLE PRECISION
  44. *          Upper bound for the smallest eigenvalue
  45. *
  46. *  SIGMA   (input/output) DOUBLE PRECISION
  47. *          Accumulated shift for the present submatrix
  48. *
  49. *  KEND    (input/output) INTEGER
  50. *          Index where minimum D(i) occurs in recurrence for
  51. *          splitting criterion
  52. *
  53. *  OFF     (input/output) INTEGER
  54. *          Offset for arrays
  55. *
  56. *  IPHASE  (input/output) INTEGER
  57. *          If IPHASE = 1 (ping) then data is in Q and E arrays
  58. *          If IPHASE = 2 (pong) then data is in QQ and EE arrays
  59. *
  60. *  ICONV   (input) INTEGER
  61. *          If ICONV = 0 a bottom part of a matrix (with a split)
  62. *          If ICONV =-3 a top part of a matrix (with a split)
  63. *
  64. *  EPS     (input) DOUBLE PRECISION
  65. *          Machine epsilon
  66. *
  67. *  TOL2    (input) DOUBLE PRECISION
  68. *          Square of the relative tolerance TOL as defined in DLASQ1
  69. *
  70. *  SMALL2  (input) DOUBLE PRECISION
  71. *          A threshold value as defined in DLASQ1
  72. *
  73. *  =====================================================================
  74. *
  75. *     .. Parameters ..
  76.       DOUBLE PRECISION   ONE, ZERO
  77.       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  78.       INTEGER            NPP
  79.       PARAMETER          ( NPP = 32 )
  80.       INTEGER            IPP
  81.       PARAMETER          ( IPP = 5 )
  82.       DOUBLE PRECISION   HALF, FOUR
  83.       PARAMETER          ( HALF = 0.5D+0, FOUR = 4.0D+0 )
  84.       INTEGER            IFLMAX
  85.       PARAMETER          ( IFLMAX = 2 )
  86. *     ..
  87. *     .. Local Scalars ..
  88.       LOGICAL            LDEF, LSPLIT
  89.       INTEGER            I, IC, ICNT, IFL, IP, ISP, K1END, K2END, KE,
  90.      $                   KS, MAXIT, N1, N2
  91.       DOUBLE PRECISION   D, DM, QEMAX, T1, TAU, TOLX, TOLY, TOLZ, XX, YY
  92. *     ..
  93. *     .. External Subroutines ..
  94.       EXTERNAL           DCOPY, DLASQ4
  95. *     ..
  96. *     .. Intrinsic Functions ..
  97.       INTRINSIC          ABS, MAX, MIN, SQRT
  98. *     ..
  99. *     .. Executable Statements ..
  100.       ICNT = 0
  101.       TAU = ZERO
  102.       DM = SUP
  103.       TOLX = SIGMA*TOL2
  104.       TOLZ = MAX( SMALL2, SIGMA )*TOL2
  105. *
  106. *     Set maximum number of iterations
  107. *
  108.       MAXIT = 100*N
  109. *
  110. *     Flipping
  111. *
  112.       IC = 2
  113.       IF( N.GT.3 ) THEN
  114.          IF( IPHASE.EQ.1 ) THEN
  115.             DO 10 I = 1, N - 2
  116.                IF( Q( I ).GT.Q( I+1 ) )
  117.      $            IC = IC + 1
  118.                IF( E( I ).GT.E( I+1 ) )
  119.      $            IC = IC + 1
  120.    10       CONTINUE
  121.             IF( Q( N-1 ).GT.Q( N ) )
  122.      $         IC = IC + 1
  123.             IF( IC.LT.N ) THEN
  124.                CALL DCOPY( N, Q, 1, QQ, -1 )
  125.                CALL DCOPY( N-1, E, 1, EE, -1 )
  126.                IF( KEND.NE.0 )
  127.      $            KEND = N - KEND + 1
  128.                IPHASE = 2
  129.             END IF
  130.          ELSE
  131.             DO 20 I = 1, N - 2
  132.                IF( QQ( I ).GT.QQ( I+1 ) )
  133.      $            IC = IC + 1
  134.                IF( EE( I ).GT.EE( I+1 ) )
  135.      $            IC = IC + 1
  136.    20       CONTINUE
  137.             IF( QQ( N-1 ).GT.QQ( N ) )
  138.      $         IC = IC + 1
  139.             IF( IC.LT.N ) THEN
  140.                CALL DCOPY( N, QQ, 1, Q, -1 )
  141.                CALL DCOPY( N-1, EE, 1, E, -1 )
  142.                IF( KEND.NE.0 )
  143.      $            KEND = N - KEND + 1
  144.                IPHASE = 1
  145.             END IF
  146.          END IF
  147.       END IF
  148.       IF( ICONV.EQ.-3 ) THEN
  149.          IF( IPHASE.EQ.1 ) THEN
  150.             GO TO 180
  151.          ELSE
  152.             GO TO 80
  153.          END IF
  154.       END IF
  155.       IF( IPHASE.EQ.2 )
  156.      $   GO TO 130
  157. *
  158. *     The ping section of the code
  159. *
  160.    30 CONTINUE
  161.       IFL = 0
  162. *
  163. *     Compute the shift
  164. *
  165.       IF( KEND.EQ.0 .OR. SUP.EQ.ZERO ) THEN
  166.          TAU = ZERO
  167.       ELSE IF( ICNT.GT.0 .AND. DM.LE.TOLZ ) THEN
  168.          TAU = ZERO
  169.       ELSE
  170.          IP = MAX( IPP, N / NPP )
  171.          N2 = 2*IP + 1
  172.          IF( N2.GE.N ) THEN
  173.             N1 = 1
  174.             N2 = N
  175.          ELSE IF( KEND+IP.GT.N ) THEN
  176.             N1 = N - 2*IP
  177.          ELSE IF( KEND-IP.LT.1 ) THEN
  178.             N1 = 1
  179.          ELSE
  180.             N1 = KEND - IP
  181.          END IF
  182.          CALL DLASQ4( N2, Q( N1 ), E( N1 ), TAU, SUP )
  183.       END IF
  184.    40 CONTINUE
  185.       ICNT = ICNT + 1
  186.       IF( ICNT.GT.MAXIT ) THEN
  187.          SUP = -ONE
  188.          RETURN
  189.       END IF
  190.       IF( TAU.EQ.ZERO ) THEN
  191. *
  192. *     dqd algorithm
  193. *
  194.          D = Q( 1 )
  195.          DM = D
  196.          KE = 0
  197.          DO 50 I = 1, N - 3
  198.             QQ( I ) = D + E( I )
  199.             D = ( D / QQ( I ) )*Q( I+1 )
  200.             IF( DM.GT.D ) THEN
  201.                DM = D
  202.                KE = I
  203.             END IF
  204.    50    CONTINUE
  205.          KE = KE + 1
  206. *
  207. *     Penultimate dqd step (in ping)
  208. *
  209.          K2END = KE
  210.          QQ( N-2 ) = D + E( N-2 )
  211.          D = ( D / QQ( N-2 ) )*Q( N-1 )
  212.          IF( DM.GT.D ) THEN
  213.             DM = D
  214.             KE = N - 1
  215.          END IF
  216. *
  217. *     Final dqd step (in ping)
  218. *
  219.          K1END = KE
  220.          QQ( N-1 ) = D + E( N-1 )
  221.          D = ( D / QQ( N-1 ) )*Q( N )
  222.          IF( DM.GT.D ) THEN
  223.             DM = D
  224.             KE = N
  225.          END IF
  226.          QQ( N ) = D
  227.       ELSE
  228. *
  229. *     The dqds algorithm (in ping)
  230. *
  231.          D = Q( 1 ) - TAU
  232.          DM = D
  233.          KE = 0
  234.          IF( D.LT.ZERO )
  235.      $      GO TO 120
  236.          DO 60 I = 1, N - 3
  237.             QQ( I ) = D + E( I )
  238.             D = ( D / QQ( I ) )*Q( I+1 ) - TAU
  239.             IF( DM.GT.D ) THEN
  240.                DM = D
  241.                KE = I
  242.                IF( D.LT.ZERO )
  243.      $            GO TO 120
  244.             END IF
  245.    60    CONTINUE
  246.          KE = KE + 1
  247. *
  248. *     Penultimate dqds step (in ping)
  249. *
  250.          K2END = KE
  251.          QQ( N-2 ) = D + E( N-2 )
  252.          D = ( D / QQ( N-2 ) )*Q( N-1 ) - TAU
  253.          IF( DM.GT.D ) THEN
  254.             DM = D
  255.             KE = N - 1
  256.             IF( D.LT.ZERO )
  257.      $         GO TO 120
  258.          END IF
  259. *
  260. *     Final dqds step (in ping)
  261. *
  262.          K1END = KE
  263.          QQ( N-1 ) = D + E( N-1 )
  264.          D = ( D / QQ( N-1 ) )*Q( N ) - TAU
  265.          IF( DM.GT.D ) THEN
  266.             DM = D
  267.             KE = N
  268.          END IF
  269.          QQ( N ) = D
  270.       END IF
  271. *
  272. *        Convergence when QQ(N) is small (in ping)
  273. *
  274.       IF( ABS( QQ( N ) ).LE.SIGMA*TOL2 ) THEN
  275.          QQ( N ) = ZERO
  276.          DM = ZERO
  277.          KE = N
  278.       END IF
  279.       IF( QQ( N ).LT.ZERO )
  280.      $   GO TO 120
  281. *
  282. *     Non-negative qd array: Update the e's
  283. *
  284.       DO 70 I = 1, N - 1
  285.          EE( I ) = ( E( I ) / QQ( I ) )*Q( I+1 )
  286.    70 CONTINUE
  287. *
  288. *     Updating sigma and iphase in ping
  289. *
  290.       SIGMA = SIGMA + TAU
  291.       IPHASE = 2
  292.    80 CONTINUE
  293.       TOLX = SIGMA*TOL2
  294.       TOLY = SIGMA*EPS
  295.       TOLZ = MAX( SIGMA, SMALL2 )*TOL2
  296. *
  297. *     Checking for deflation and convergence (in ping)
  298. *
  299.    90 CONTINUE
  300.       IF( N.LE.2 )
  301.      $   RETURN
  302. *
  303. *        Deflation: bottom 1x1 (in ping)
  304. *
  305.       LDEF = .FALSE.
  306.       IF( EE( N-1 ).LE.TOLZ ) THEN
  307.          LDEF = .TRUE.
  308.       ELSE IF( SIGMA.GT.ZERO ) THEN
  309.          IF( EE( N-1 ).LE.EPS*( SIGMA+QQ( N ) ) ) THEN
  310.             IF( EE( N-1 )*( QQ( N ) / ( QQ( N )+SIGMA ) ).LE.TOL2*
  311.      $          ( QQ( N )+SIGMA ) ) THEN
  312.                LDEF = .TRUE.
  313.             END IF
  314.          END IF
  315.       ELSE
  316.          IF( EE( N-1 ).LE.QQ( N )*TOL2 ) THEN
  317.             LDEF = .TRUE.
  318.          END IF
  319.       END IF
  320.       IF( LDEF ) THEN
  321.          Q( N ) = QQ( N ) + SIGMA
  322.          N = N - 1
  323.          ICONV = ICONV + 1
  324.          GO TO 90
  325.       END IF
  326. *
  327. *        Deflation: bottom 2x2 (in ping)
  328. *
  329.       LDEF = .FALSE.
  330.       IF( EE( N-2 ).LE.TOLZ ) THEN
  331.          LDEF = .TRUE.
  332.       ELSE IF( SIGMA.GT.ZERO ) THEN
  333.          T1 = SIGMA + EE( N-1 )*( SIGMA / ( SIGMA+QQ( N ) ) )
  334.          IF( EE( N-2 )*( T1 / ( QQ( N-1 )+T1 ) ).LE.TOLY ) THEN
  335.             IF( EE( N-2 )*( QQ( N-1 ) / ( QQ( N-1 )+T1 ) ).LE.TOLX )
  336.      $           THEN
  337.                LDEF = .TRUE.
  338.             END IF
  339.          END IF
  340.       ELSE
  341.          IF( EE( N-2 ).LE.( QQ( N ) / ( QQ( N )+EE( N-1 )+QQ( N-1 ) ) )*
  342.      $       QQ( N-1 )*TOL2 ) THEN
  343.             LDEF = .TRUE.
  344.          END IF
  345.       END IF
  346.       IF( LDEF ) THEN
  347.          QEMAX = MAX( QQ( N ), QQ( N-1 ), EE( N-1 ) )
  348.          IF( QEMAX.NE.ZERO ) THEN
  349.             IF( QEMAX.EQ.QQ( N-1 ) ) THEN
  350.                XX = HALF*( QQ( N )+QQ( N-1 )+EE( N-1 )+QEMAX*
  351.      $              SQRT( ( ( QQ( N )-QQ( N-1 )+EE( N-1 ) ) /
  352.      $              QEMAX )**2+FOUR*EE( N-1 ) / QEMAX ) )
  353.             ELSE IF( QEMAX.EQ.QQ( N ) ) THEN
  354.                XX = HALF*( QQ( N )+QQ( N-1 )+EE( N-1 )+QEMAX*
  355.      $              SQRT( ( ( QQ( N-1 )-QQ( N )+EE( N-1 ) ) /
  356.      $              QEMAX )**2+FOUR*EE( N-1 ) / QEMAX ) )
  357.             ELSE
  358.                XX = HALF*( QQ( N )+QQ( N-1 )+EE( N-1 )+QEMAX*
  359.      $              SQRT( ( ( QQ( N )-QQ( N-1 )+EE( N-1 ) ) /
  360.      $              QEMAX )**2+FOUR*QQ( N-1 ) / QEMAX ) )
  361.             END IF
  362.             YY = ( MAX( QQ( N ), QQ( N-1 ) ) / XX )*
  363.      $           MIN( QQ( N ), QQ( N-1 ) )
  364.          ELSE
  365.             XX = ZERO
  366.             YY = ZERO
  367.          END IF
  368.          Q( N-1 ) = SIGMA + XX
  369.          Q( N ) = YY + SIGMA
  370.          N = N - 2
  371.          ICONV = ICONV + 2
  372.          GO TO 90
  373.       END IF
  374. *
  375. *     Updating bounds before going to pong
  376. *
  377.       IF( ICONV.EQ.0 ) THEN
  378.          KEND = KE
  379.          SUP = MIN( DM, SUP-TAU )
  380.       ELSE IF( ICONV.GT.0 ) THEN
  381.          SUP = MIN( QQ( N ), QQ( N-1 ), QQ( N-2 ), QQ( 1 ), QQ( 2 ),
  382.      $         QQ( 3 ) )
  383.          IF( ICONV.EQ.1 ) THEN
  384.             KEND = K1END
  385.          ELSE IF( ICONV.EQ.2 ) THEN
  386.             KEND = K2END
  387.          ELSE
  388.             KEND = N
  389.          END IF
  390.          ICNT = 0
  391.          MAXIT = 100*N
  392.       END IF
  393. *
  394. *     Checking for splitting in ping
  395. *
  396.       LSPLIT = .FALSE.
  397.       DO 100 KS = N - 3, 3, -1
  398.          IF( EE( KS ).LE.TOLY ) THEN
  399.             IF( EE( KS )*( MIN( QQ( KS+1 ),
  400.      $          QQ( KS ) ) / ( MIN( QQ( KS+1 ), QQ( KS ) )+SIGMA ) ).LE.
  401.      $          TOLX ) THEN
  402.                LSPLIT = .TRUE.
  403.                GO TO 110
  404.             END IF
  405.          END IF
  406.   100 CONTINUE
  407. *
  408.       KS = 2
  409.       IF( EE( 2 ).LE.TOLZ ) THEN
  410.          LSPLIT = .TRUE.
  411.       ELSE IF( SIGMA.GT.ZERO ) THEN
  412.          T1 = SIGMA + EE( 1 )*( SIGMA / ( SIGMA+QQ( 1 ) ) )
  413.          IF( EE( 2 )*( T1 / ( QQ( 1 )+T1 ) ).LE.TOLY ) THEN
  414.             IF( EE( 2 )*( QQ( 1 ) / ( QQ( 1 )+T1 ) ).LE.TOLX ) THEN
  415.                LSPLIT = .TRUE.
  416.             END IF
  417.          END IF
  418.       ELSE
  419.          IF( EE( 2 ).LE.( QQ( 1 ) / ( QQ( 1 )+EE( 1 )+QQ( 2 ) ) )*
  420.      $       QQ( 2 )*TOL2 ) THEN
  421.             LSPLIT = .TRUE.
  422.          END IF
  423.       END IF
  424.       IF( LSPLIT )
  425.      $   GO TO 110
  426. *
  427.       KS = 1
  428.       IF( EE( 1 ).LE.TOLZ ) THEN
  429.          LSPLIT = .TRUE.
  430.       ELSE IF( SIGMA.GT.ZERO ) THEN
  431.          IF( EE( 1 ).LE.EPS*( SIGMA+QQ( 1 ) ) ) THEN
  432.             IF( EE( 1 )*( QQ( 1 ) / ( QQ( 1 )+SIGMA ) ).LE.TOL2*
  433.      $          ( QQ( 1 )+SIGMA ) ) THEN
  434.                LSPLIT = .TRUE.
  435.             END IF
  436.          END IF
  437.       ELSE
  438.          IF( EE( 1 ).LE.QQ( 1 )*TOL2 ) THEN
  439.             LSPLIT = .TRUE.
  440.          END IF
  441.       END IF
  442. *
  443.   110 CONTINUE
  444.       IF( LSPLIT ) THEN
  445.          SUP = MIN( QQ( N ), QQ( N-1 ), QQ( N-2 ) )
  446.          ISP = -( OFF+1 )
  447.          OFF = OFF + KS
  448.          N = N - KS
  449.          KEND = MAX( 1, KEND-KS )
  450.          E( KS ) = SIGMA
  451.          EE( KS ) = ISP
  452.          ICONV = 0
  453.          RETURN
  454.       END IF
  455. *
  456. *     Coincidence
  457. *
  458.       IF( TAU.EQ.ZERO .AND. DM.LE.TOLZ .AND. KEND.NE.N .AND. ICONV.EQ.
  459.      $    0 .AND. ICNT.GT.0 ) THEN
  460.          CALL DCOPY( N-KE, E( KE ), 1, QQ( KE ), 1 )
  461.          QQ( N ) = ZERO
  462.          CALL DCOPY( N-KE, Q( KE+1 ), 1, EE( KE ), 1 )
  463.          SUP = ZERO
  464.       END IF
  465.       ICONV = 0
  466.       GO TO 130
  467. *
  468. *     A new shift when the previous failed (in ping)
  469. *
  470.   120 CONTINUE
  471.       IFL = IFL + 1
  472.       SUP = TAU
  473. *
  474. *     SUP is small or
  475. *     Too many bad shifts (ping)
  476. *
  477.       IF( SUP.LE.TOLZ .OR. IFL.GE.IFLMAX ) THEN
  478.          TAU = ZERO
  479.          GO TO 40
  480. *
  481. *     The asymptotic shift (in ping)
  482. *
  483.       ELSE
  484.          TAU = MAX( TAU+D, ZERO )
  485.          IF( TAU.LE.TOLZ )
  486.      $      TAU = ZERO
  487.          GO TO 40
  488.       END IF
  489. *
  490. *     the pong section of the code
  491. *
  492.   130 CONTINUE
  493.       IFL = 0
  494. *
  495. *     Compute the shift (in pong)
  496. *
  497.       IF( KEND.EQ.0 .AND. SUP.EQ.ZERO ) THEN
  498.          TAU = ZERO
  499.       ELSE IF( ICNT.GT.0 .AND. DM.LE.TOLZ ) THEN
  500.          TAU = ZERO
  501.       ELSE
  502.          IP = MAX( IPP, N / NPP )
  503.          N2 = 2*IP + 1
  504.          IF( N2.GE.N ) THEN
  505.             N1 = 1
  506.             N2 = N
  507.          ELSE IF( KEND+IP.GT.N ) THEN
  508.             N1 = N - 2*IP
  509.          ELSE IF( KEND-IP.LT.1 ) THEN
  510.             N1 = 1
  511.          ELSE
  512.             N1 = KEND - IP
  513.          END IF
  514.          CALL DLASQ4( N2, QQ( N1 ), EE( N1 ), TAU, SUP )
  515.       END IF
  516.   140 CONTINUE
  517.       ICNT = ICNT + 1
  518.       IF( ICNT.GT.MAXIT ) THEN
  519.          SUP = -SUP
  520.          RETURN
  521.       END IF
  522.       IF( TAU.EQ.ZERO ) THEN
  523. *
  524. *     The dqd algorithm (in pong)
  525. *
  526.          D = QQ( 1 )
  527.          DM = D
  528.          KE = 0
  529.          DO 150 I = 1, N - 3
  530.             Q( I ) = D + EE( I )
  531.             D = ( D / Q( I ) )*QQ( I+1 )
  532.             IF( DM.GT.D ) THEN
  533.                DM = D
  534.                KE = I
  535.             END IF
  536.   150    CONTINUE
  537.          KE = KE + 1
  538. *
  539. *     Penultimate dqd step (in pong)
  540. *
  541.          K2END = KE
  542.          Q( N-2 ) = D + EE( N-2 )
  543.          D = ( D / Q( N-2 ) )*QQ( N-1 )
  544.          IF( DM.GT.D ) THEN
  545.             DM = D
  546.             KE = N - 1
  547.          END IF
  548. *
  549. *     Final dqd step (in pong)
  550. *
  551.          K1END = KE
  552.          Q( N-1 ) = D + EE( N-1 )
  553.          D = ( D / Q( N-1 ) )*QQ( N )
  554.          IF( DM.GT.D ) THEN
  555.             DM = D
  556.             KE = N
  557.          END IF
  558.          Q( N ) = D
  559.       ELSE
  560. *
  561. *     The dqds algorithm (in pong)
  562. *
  563.          D = QQ( 1 ) - TAU
  564.          DM = D
  565.          KE = 0
  566.          IF( D.LT.ZERO )
  567.      $      GO TO 220
  568.          DO 160 I = 1, N - 3
  569.             Q( I ) = D + EE( I )
  570.             D = ( D / Q( I ) )*QQ( I+1 ) - TAU
  571.             IF( DM.GT.D ) THEN
  572.                DM = D
  573.                KE = I
  574.                IF( D.LT.ZERO )
  575.      $            GO TO 220
  576.             END IF
  577.   160    CONTINUE
  578.          KE = KE + 1
  579. *
  580. *     Penultimate dqds step (in pong)
  581. *
  582.          K2END = KE
  583.          Q( N-2 ) = D + EE( N-2 )
  584.          D = ( D / Q( N-2 ) )*QQ( N-1 ) - TAU
  585.          IF( DM.GT.D ) THEN
  586.             DM = D
  587.             KE = N - 1
  588.             IF( D.LT.ZERO )
  589.      $         GO TO 220
  590.          END IF
  591. *
  592. *     Final dqds step (in pong)
  593. *
  594.          K1END = KE
  595.          Q( N-1 ) = D + EE( N-1 )
  596.          D = ( D / Q( N-1 ) )*QQ( N ) - TAU
  597.          IF( DM.GT.D ) THEN
  598.             DM = D
  599.             KE = N
  600.          END IF
  601.          Q( N ) = D
  602.       END IF
  603. *
  604. *        Convergence when is small (in pong)
  605. *
  606.       IF( ABS( Q( N ) ).LE.SIGMA*TOL2 ) THEN
  607.          Q( N ) = ZERO
  608.          DM = ZERO
  609.          KE = N
  610.       END IF
  611.       IF( Q( N ).LT.ZERO )
  612.      $   GO TO 220
  613. *
  614. *     Non-negative qd array: Update the e's
  615. *
  616.       DO 170 I = 1, N - 1
  617.          E( I ) = ( EE( I ) / Q( I ) )*QQ( I+1 )
  618.   170 CONTINUE
  619. *
  620. *     Updating sigma and iphase in pong
  621. *
  622.       SIGMA = SIGMA + TAU
  623.   180 CONTINUE
  624.       IPHASE = 1
  625.       TOLX = SIGMA*TOL2
  626.       TOLY = SIGMA*EPS
  627. *
  628. *     Checking for deflation and convergence (in pong)
  629. *
  630.   190 CONTINUE
  631.       IF( N.LE.2 )
  632.      $   RETURN
  633. *
  634. *        Deflation: bottom 1x1 (in pong)
  635. *
  636.       LDEF = .FALSE.
  637.       IF( E( N-1 ).LE.TOLZ ) THEN
  638.          LDEF = .TRUE.
  639.       ELSE IF( SIGMA.GT.ZERO ) THEN
  640.          IF( E( N-1 ).LE.EPS*( SIGMA+Q( N ) ) ) THEN
  641.             IF( E( N-1 )*( Q( N ) / ( Q( N )+SIGMA ) ).LE.TOL2*
  642.      $          ( Q( N )+SIGMA ) ) THEN
  643.                LDEF = .TRUE.
  644.             END IF
  645.          END IF
  646.       ELSE
  647.          IF( E( N-1 ).LE.Q( N )*TOL2 ) THEN
  648.             LDEF = .TRUE.
  649.          END IF
  650.       END IF
  651.       IF( LDEF ) THEN
  652.          Q( N ) = Q( N ) + SIGMA
  653.          N = N - 1
  654.          ICONV = ICONV + 1
  655.          GO TO 190
  656.       END IF
  657. *
  658. *        Deflation: bottom 2x2 (in pong)
  659. *
  660.       LDEF = .FALSE.
  661.       IF( E( N-2 ).LE.TOLZ ) THEN
  662.          LDEF = .TRUE.
  663.       ELSE IF( SIGMA.GT.ZERO ) THEN
  664.          T1 = SIGMA + E( N-1 )*( SIGMA / ( SIGMA+Q( N ) ) )
  665.          IF( E( N-2 )*( T1 / ( Q( N-1 )+T1 ) ).LE.TOLY ) THEN
  666.             IF( E( N-2 )*( Q( N-1 ) / ( Q( N-1 )+T1 ) ).LE.TOLX ) THEN
  667.                LDEF = .TRUE.
  668.             END IF
  669.          END IF
  670.       ELSE
  671.          IF( E( N-2 ).LE.( Q( N ) / ( Q( N )+EE( N-1 )+Q( N-1 ) )*Q( N-
  672.      $       1 ) )*TOL2 ) THEN
  673.             LDEF = .TRUE.
  674.          END IF
  675.       END IF
  676.       IF( LDEF ) THEN
  677.          QEMAX = MAX( Q( N ), Q( N-1 ), E( N-1 ) )
  678.          IF( QEMAX.NE.ZERO ) THEN
  679.             IF( QEMAX.EQ.Q( N-1 ) ) THEN
  680.                XX = HALF*( Q( N )+Q( N-1 )+E( N-1 )+QEMAX*
  681.      $              SQRT( ( ( Q( N )-Q( N-1 )+E( N-1 ) ) / QEMAX )**2+
  682.      $              FOUR*E( N-1 ) / QEMAX ) )
  683.             ELSE IF( QEMAX.EQ.Q( N ) ) THEN
  684.                XX = HALF*( Q( N )+Q( N-1 )+E( N-1 )+QEMAX*
  685.      $              SQRT( ( ( Q( N-1 )-Q( N )+E( N-1 ) ) / QEMAX )**2+
  686.      $              FOUR*E( N-1 ) / QEMAX ) )
  687.             ELSE
  688.                XX = HALF*( Q( N )+Q( N-1 )+E( N-1 )+QEMAX*
  689.      $              SQRT( ( ( Q( N )-Q( N-1 )+E( N-1 ) ) / QEMAX )**2+
  690.      $              FOUR*Q( N-1 ) / QEMAX ) )
  691.             END IF
  692.             YY = ( MAX( Q( N ), Q( N-1 ) ) / XX )*
  693.      $           MIN( Q( N ), Q( N-1 ) )
  694.          ELSE
  695.             XX = ZERO
  696.             YY = ZERO
  697.          END IF
  698.          Q( N-1 ) = SIGMA + XX
  699.          Q( N ) = YY + SIGMA
  700.          N = N - 2
  701.          ICONV = ICONV + 2
  702.          GO TO 190
  703.       END IF
  704. *
  705. *     Updating bounds before going to pong
  706. *
  707.       IF( ICONV.EQ.0 ) THEN
  708.          KEND = KE
  709.          SUP = MIN( DM, SUP-TAU )
  710.       ELSE IF( ICONV.GT.0 ) THEN
  711.          SUP = MIN( Q( N ), Q( N-1 ), Q( N-2 ), Q( 1 ), Q( 2 ), Q( 3 ) )
  712.          IF( ICONV.EQ.1 ) THEN
  713.             KEND = K1END
  714.          ELSE IF( ICONV.EQ.2 ) THEN
  715.             KEND = K2END
  716.          ELSE
  717.             KEND = N
  718.          END IF
  719.          ICNT = 0
  720.          MAXIT = 100*N
  721.       END IF
  722. *
  723. *     Checking for splitting in pong
  724. *
  725.       LSPLIT = .FALSE.
  726.       DO 200 KS = N - 3, 3, -1
  727.          IF( E( KS ).LE.TOLY ) THEN
  728.             IF( E( KS )*( MIN( Q( KS+1 ), Q( KS ) ) / ( MIN( Q( KS+1 ),
  729.      $          Q( KS ) )+SIGMA ) ).LE.TOLX ) THEN
  730.                LSPLIT = .TRUE.
  731.                GO TO 210
  732.             END IF
  733.          END IF
  734.   200 CONTINUE
  735. *
  736.       KS = 2
  737.       IF( E( 2 ).LE.TOLZ ) THEN
  738.          LSPLIT = .TRUE.
  739.       ELSE IF( SIGMA.GT.ZERO ) THEN
  740.          T1 = SIGMA + E( 1 )*( SIGMA / ( SIGMA+Q( 1 ) ) )
  741.          IF( E( 2 )*( T1 / ( Q( 1 )+T1 ) ).LE.TOLY ) THEN
  742.             IF( E( 2 )*( Q( 1 ) / ( Q( 1 )+T1 ) ).LE.TOLX ) THEN
  743.                LSPLIT = .TRUE.
  744.             END IF
  745.          END IF
  746.       ELSE
  747.          IF( E( 2 ).LE.( Q( 1 ) / ( Q( 1 )+E( 1 )+Q( 2 ) ) )*Q( 2 )*
  748.      $       TOL2 ) THEN
  749.             LSPLIT = .TRUE.
  750.          END IF
  751.       END IF
  752.       IF( LSPLIT )
  753.      $   GO TO 210
  754. *
  755.       KS = 1
  756.       IF( E( 1 ).LE.TOLZ ) THEN
  757.          LSPLIT = .TRUE.
  758.       ELSE IF( SIGMA.GT.ZERO ) THEN
  759.          IF( E( 1 ).LE.EPS*( SIGMA+Q( 1 ) ) ) THEN
  760.             IF( E( 1 )*( Q( 1 ) / ( Q( 1 )+SIGMA ) ).LE.TOL2*
  761.      $          ( Q( 1 )+SIGMA ) ) THEN
  762.                LSPLIT = .TRUE.
  763.             END IF
  764.          END IF
  765.       ELSE
  766.          IF( E( 1 ).LE.Q( 1 )*TOL2 ) THEN
  767.             LSPLIT = .TRUE.
  768.          END IF
  769.       END IF
  770. *
  771.   210 CONTINUE
  772.       IF( LSPLIT ) THEN
  773.          SUP = MIN( Q( N ), Q( N-1 ), Q( N-2 ) )
  774.          ISP = OFF + 1
  775.          OFF = OFF + KS
  776.          KEND = MAX( 1, KEND-KS )
  777.          N = N - KS
  778.          E( KS ) = SIGMA
  779.          EE( KS ) = ISP
  780.          ICONV = 0
  781.          RETURN
  782.       END IF
  783. *
  784. *     Coincidence
  785. *
  786.       IF( TAU.EQ.ZERO .AND. DM.LE.TOLZ .AND. KEND.NE.N .AND. ICONV.EQ.
  787.      $    0 .AND. ICNT.GT.0 ) THEN
  788.          CALL DCOPY( N-KE, EE( KE ), 1, Q( KE ), 1 )
  789.          Q( N ) = ZERO
  790.          CALL DCOPY( N-KE, QQ( KE+1 ), 1, E( KE ), 1 )
  791.          SUP = ZERO
  792.       END IF
  793.       ICONV = 0
  794.       GO TO 30
  795. *
  796. *     Computation of a new shift when the previous failed (in pong)
  797. *
  798.   220 CONTINUE
  799.       IFL = IFL + 1
  800.       SUP = TAU
  801. *
  802. *     SUP is small or
  803. *     Too many bad shifts (in pong)
  804. *
  805.       IF( SUP.LE.TOLZ .OR. IFL.GE.IFLMAX ) THEN
  806.          TAU = ZERO
  807.          GO TO 140
  808. *
  809. *     The asymptotic shift (in pong)
  810. *
  811.       ELSE
  812.          TAU = MAX( TAU+D, ZERO )
  813.          IF( TAU.LE.TOLZ )
  814.      $      TAU = ZERO
  815.          GO TO 140
  816.       END IF
  817. *
  818. *     End of DLASQ3
  819. *
  820.       END
  821.